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Abstract. We present a new approach for statistical inference on noise properties of CMB anisotropy data. 
We consider a Maximum Likelihood parametric estimator to recover the full dependence structure of the noise 
process. We also consider a semiparametric procedure which is only sensitive to the low frequency behavior of 
the noise spectral density. Both approaches are statistically robust and computationally convenient in the case of 
long memory noise, even under nonstationary circumstances. We show that noise properties can be consistently 
derived by such procedures without resorting to currently used iterative noise-signal methods. More importantly, 
we show that optimal (GLS) CMB maps can be obtained from the observed timestream with the only knowledge 
of the noise memory parameter, the outcome of our estimators. 
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1. Introduction 

Recent high sensitivity, high resolution observations of 
the Cosmic Microwave Background (CMB) anisotropy are 
posing strong constraints on Cosmology (see, e.g., Balbi 
et al., 2000; de Bernardis et al., 2001). The ability to dis- 
criminate among different cosmological scenarios critically 
relies on the efficiency and accuracy of the data analysis 
pipeline, other than -obviously- on experimental know- 
how. In order to be efficient, data analysis has to deal with 
the size of the data sets, which is bound to increase with 
new space missions such as MAPQ and Planck Surveyor^. 
In order to be accurate, the data analysis pipeline should 
implement optimal statistical techniques to reconstruct 
the correlation structure of receiver noise, i.e. the behavior 
of its spectral density. This is a critical point for data anal- 
ysis of CMB anisotropy one- horned experiment. In fact, 
both Maximum Likelihood map making (see e.g. Wright, 
1996; Borrill, 1999; Natofi et al., 2001) and angular power 
spectrum estimation (Tegmark 1997; Bond et al., 1998; 
Oh et al., 1999) heavily depend on the knowledge of the 
noise spectral density, which of course must be determined 
from the data themselves. 



Our paper is not the first to focus on estimating the 
noise properties out of a CMB anisotropy data stream. 
Ferreira & Jaffe (2000) advocated simultaneous estimation 
of signal and noise on the basis of a Bayesian approach. 
However, their estimates for the noise spectral density are 
piecewise constant with several discontinuities. Far from 
being realistic, this assumption does not apply to long 
range dependence models. Prunet et al. (2001) obtained a 
purely nonparametric estimate of the noise spectral den- 
sity as a by-product of their iterative map making pro- 
cedure. Even if cast in a frequentist framework, the re- 
sulting algorithm is very similar to Ferreira and Jaffe's 
proposal. It is worth noting that, in both cases, purely 
nonparametric estimates of noise properties lack a rigor- 
ous justification in the statistical literature if noise is a 
long memory process. In this case, all of these difficulties 
can be overcome by using well known and robust statisti- 
cal techniques which are well supported by the most recent 
literature (Fox & Taqqu, 1986; Dahlhaus, 1995; Robinson 
1995ab; Giraitis & Taqqu, 1999). In Sect. 2 we discuss 
these techniques in the context of CMB data analysis. 
In Sect. 3 we apply them to simulated datasets of the 
Planck and BOOOMERanG| experiments. Finally, in 
Sect. 4 we draw our main conclusions and point out direc- 
tions for future research. 
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2. EfFicient estimation of tlie noise properties 

We assume that the data, d, consist of the sum of two 
uncorrelated time series: 



St + nt 



(1) 



for f = 1, 2, . . . , T. Here n and s are noise and signal, re- 
spectively. The latter is modeled assuming linear depen- 
dence on the sky pattern: 



St 



K 

PtpTOp 



(2) 



We imagine the "map" m to be smeared by the instru- 
mental beam and discretized into K elements (or pixels). 
The matrix P effectively defines the experimental scan- 
ning strategy by "unrolling" the sky map into a time se- 
ries. We limit ourselves to the case of one-horned experi- 
ments and symmetric beam profile, for which the matrix 
P has a very simple structure, consisting of a single non 
zero entry (with value 1) per row. This entry corresponds 
to the pixel observed at a given time t. For very low signal 
to noise (S/N) ratio the timestream, d, is noise dominated. 
In this Section we pretend that this is the case. We will 
show below, in Sect. 3, how to drop this assumption. 

2.1. Long memory processes 

A process is said to be wide sense (or second order) sta- 
tionary if its mean, variance and covariance functions are 
finite and constant over time: {nt) — /i, (ntUt) — 7(0) < 

00 and (ntnt-T) = i{t), t — ±1,±2, For a noise 

process we shall take /i = 0. The Cramer Representation 
Theorem ensures that for any second order stationary pro- 
cess we can write: 



nt = 



eM^tf)VPif)dWif) 



(3) 



where W{f) is a random measure such that {dW{f)) — 0, 
while (dW {f)dW (g)) = df \i f = g and otherwise (i.e. 
W represents a white noise measure). The function P{f) 
is the spectral density function of the process n, uniquely 
related to the sequence of covariances 7(t) through the 
usual Fourier formulae: 



ZTT — ^ 

T — — 00 

7(r) = r exp{-tTf)P{f )df . 

J — TT 

Thus, wide sense stationarity implies: 

7(0) = / P{f)df < 00 , 



that is, the spectral density needs always be integrable. 
State of the art microwave detectors produce "1//" noise. 
We shall therefore focus on long memory (or long range 



dependent) noise. For this processes the spectral density 
is such that 



lim^ 



(4) 



where < G < 00 and < a < 1. So, even if divergent at 
zero frequency, the spectral density is always integrable. 
Eq. (0) can be discretized as follows: 



1 



nt 



(5) 



Here fj — 2Trj/T are Fourier frequencies, while ej denotes 
a sequence of zero-mean, independent and identically dis- 
tributed (i.i.d.) variables with unit variance. Eq. (^ is 
widely used to generate random processes with a given 
spectral density P{f) in Monte Carlo procedures. If we 
take, for instance, P{f) — Gf^", straightforward algebra 
shows that: 



1 

lim (ntnt) ~ lim __ ^P(/j) = 

.7 = 1 



0(1), a < 1 
O(logr), a=l 
0{T"~^), a > 1 



Therefore, if a > 1 the variance of n grows with time, i.e. 
the process is not stationary. We shall first focus on the 
stationary case a < 1. This assumption will be relaxed to 



a < 2 in Sect. 2.4 



2.2. Fully parametric procedures 

Assume first that the the functional form of the spectral 
density for the process n is known and completely specified 
by a set of parameters 9 = (6'i, 6*2, . . . , 6*^). For "1//" 
noise it is customary to write: 



P{f, 



Ail 



(6) 



9 = {A,a,fk), where A is an amplitude, fk is the so- 
called knee frequency and a is the spectral index intro- 
duced above. Throughout this paper we shall assume that 
n is Gaussian even if this condition is not essential for our 
conclusions to hold. Under Gaussianity n has a multivari- 
ate density function: 



*(n) = 



exp 



-n'N^^n 



(7) 



where N is the noise covariance matrix and the prime (') 
denotes transposition. Now, for large T the matrix N^^ 
is nearly circulant and can be approximated as: 



(8) 



where A = diag{P{fi), ...P{fT/2)}^ the dagger (f ) denotes 
complex conjugation coupled with transposition and Q is 
the rectangular (T/2 x T) matrix 



(9) 
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The sample Discrete Fourier Transform (DFT) 



2.3. Semiparametric procedures 



•n = {wTifi), ■ ■ • ,wt(/t/2)} , 



has elements 



The sample periodogram is defined as 

lT{fi)^WT{f,)wT{f,)* , 

where the asterisk denotes conjugation. 

From Eq. (|^), it is immediate to derive the log- 
likelihood function for as 

C{e- n) = log2^ - i log |N(0)| - in'N(0)-in . 

Maximization of this form in the time domain is a compu- 
tationally unfeasible task even for moderately large values 
of T, let aside T ^ 10^ that will become customary in 
CMB experiments. Based only on the approximations in 
Eq. (|^) and building upon ideas that trace back to work 
by Whittle (1953), Fox and Taqqu (1986) have suggested 
the following frequency domain expression, valid up to a 
constant term: 

-2L{e-n) ~ log|QAQt|+n'QtA-iQn 

= log|QtQ|+iog|A| + (Qn)U-iQn 

T/2 



(10) 



The Whittle Estimate (WE) for the noise parameters, 6, 
is then defined as: 



T/2 



T/2 



e = arg nrin I - ^ log P(/, ; 0) - ^ 

where argmina;{/(x)} denotes the value x which mini- 
mizes /(•). The same authors show that, under regularity 
conditions, these estimates are Vr-consistent, asymptot- 
ically unbiased and Gaussian, meaning that 

lim \/r(^-6l) -A/'(0,S]) (11) 

T — >OQ 

where J\f denotes a Gaussian multivariate with covariance 
matrix S whose explicit expression is given by Fox and 
Taqqu (1986). The estimate 9 is computationally very 
convenient and Dahlhaus (1989) proves that it is asymp- 
totically equivalent to the corresponding estimate in the 
time domain, i.e. to standard Maximum Likelihood. It fol- 
lows immediately that the WE is absolutely efficient in the 
Cramer-Rao sense, that is it achieves minimum variance. 
Generalizations to a non Gaussian framework are con- 
sidered by Giraitis and Surgailis (1990) and Giraitis and 
Taqqu (1999). Thus, provided that the functional form 
of the spectral density is known a priori, WE solves the 
problem of optimal inference of the noise properties given 
that the observed timestream is noise dominated. 



In some situations the functional form of the spectral den- 
sity may not be known. Nonetheless, in the presence of 
long range dependence, knowledge of the low frequency be- 
havior of the noise spectral density is all that is needed to 
implement optimal filtering procedures (Dahlhaus, 1995). 
Thus, determining a is sufficient for Generalized Least 
Squares (GLS) or Maximum Likelihood (ML) map mak- 
ing. This fact calls for statistical methods which, rather 
than parameterize the full spectral density, only rely on 
the much milder condition of Eq. (^). One such procedure, 
the so-called Log Periodogram Estimate (LPE), was intro- 
duced by Geweke and Porter-Hudak (1983) and discussed 
rigorously by Robinson (1995a). Consider the identity: 



log It (/, ) = log G - a log fj -f log 



log 



ItUjI 
(12) 



where j — l,2,...,m < T/2 and m = m{T) is chosen 
such that limy^oo 1/m = and limT_,oo rn/T = 0. Under 
this Asymptotic Bandwidth Condition (ABC) and in view 



of Eq. 



one has that log[P(/j)/G/: 



whereas 



\og[lTifj)/Pifj)] can be approximated as a sequence of 
zero-mean and nearly uncorrelated residuals, with finite 
variance. Therefore, it seems natural to consider Eq. (12) 
as a regression model. This heuristic is made rigorous 
by Robinson (1995a), where the Ordinary Least Square 
(OLS) estimate is considered: 



argmin V {log It (/j) 



G-«log/,}= 



(13) 



The intuitive meaning of the ABC is that we only consider 
a vanishing (as T — *■ oo) subset of frequencies around the 
origin. Because of this, we are actually discarding most 
of the available information as it is the case for semipara- 
metric procedures. The practical consequence, highlighted 
in Robinson (1995a), is that a, despite being asymptoti- 
cally Gaussian and unbiased, is only y^-consistent (in the 
sense of Eq. ([III)). It has therefore asymptotical efficiency 
zero with respect to WE when the spectral density is cor- 
rectly parameterized. In principle, a misspecified model 
may entail inconsistent estimates of all parameters, and 
in particular of a, which is the parameter of interest for 
many, if not most, applications. In this context, robust 
semiparametric procedures may represent a more reliable 
choice. 

Robinson (1995b) considers another semiparametric 
method which can be viewed as a narrow-band analogous 
of WE. Its properties and motivating rationale, however, 
are too close to LPE to warrant independent consideration 
in this paper. 

A somewhat intermediate attitude between the fully 
parametric WE and the semiparametric LPE has been 
recently set forth by Moulines & Soulier (1999). The idea 
is basically to use Eq. (^ to analyze the lowest frequencies, 
where the long range dependent behavior shows up, and to 
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consider a series expansion into orthogonal components for 
the remaining part of the spectral density. These estimates 
are in principle appealing because they can be shown to 
converge as fast as \/T j logT, provided that the spectral 
density is sufficiently well behaved close to the origin 
(see Moulines & Soulier, 1999) for more details). However, 
they are computationally more costly than WE and LPE 
because they require multiple regression with hundreds or 
thousands of estimated explanatory terms. 

2 A. The nonstationary case 

The analysis in the previous Section assumes stationarity 
of the noise series. This assumption may turn out to be too 
strong over a long time span. For instance, the variance of 
the noise might grow (slowly) as the observation goes on. 
A general model for such a process might be 

t 

where St represents a sequence of independent innovations 
with zero mean and finite variance, whereas the weight 
sequence ~ kj"/-^"^, (k > 0). For a < 1, it is not 
difficult to see that n is asymptotically stationary and 
satisfies Eq. (|) as t ^ oo. For 1 < a < 2 Velasco (1999) 
proves that, although the spectral density is no longer de- 
fined, the expected value of the periodogram maintains 
the same behavior as in Eq. (Q). In the case a = 2, n 
behaves (asymptotically) as a pure random walk process: 

t 

j=o 

The case a > 2 can also be considered. Robinson & 
Marinucci (2001) investigate the behavior of the peri- 
odogram for any positive value of a. However, this would 
entail a superlinear growth of the variance with time. We 
rule this out as experimentally questionable. 

LPE and WE are shown to yield consistent estimates 
even under the condition a < 2 (Velasco, 1999; Robinson 
& Velasco, 2000). This fact allows us to also consider val- 
ues of a greater than unity in the simulations we will carry 
out in the next Section. As a consequence, the procedures 
advocated in this paper enjoy a marked advantage with re- 
spect to those so far considered in the CMB literature, to 
the extent in which the latter analyses are deeply rooted 
in a stationary framework and apparently lack rigorous 
justification otherwise. 

3. Numerical Applications 

The methods outlined above assume that the noise time 
series, n, is known. Obviously this is not the case for a real 
world experiment where the observed sample is a combi- 
nation of noise and signal. This is precisely the reason why 
iterative methods have grown attention in the literature. 



In this paper, however, we will follow a different - yet sim- 
pler - approach. That is, we will test the efficiency of the 
statistical methods outlined above on the simplest noise 
estimator, obtained by subtracting the OLS (i.e. naively 
coadded) map from the data. To be more specific, let us 
write Eq. (Q) in matrix form as: 

d = Pm -I- n. 

The OLS estimator is then: 

n = d - P(P^P)"ip^d = [I - P(P^P)-ip'^]n (14) 

which does not include any contribution from the signal, 
independently of the specific S/N ratio of the timestrcam 
(Natoh et al. 2001). 

The plan of this Section is as follows. First, we test 
the efficiency of WE and LPE on the pure noise (n) 
timestrcam. We then compare these results with those 
obtained by using the n estimator. To assess this point, 
we make Monte Carlo simulations for the Planck mis- 
sion. We then use our noise estimates to perform map 
making on simulated datasets for both Planck and 
BOOMERanG. In doing so, we also compare the outcome 
of our noise estimation method with our implementation 
of the iterative scheme proposed by Prunet et al. (2001) 
for BOOMERanG. 

3.1. LPE and WE efficiency on n and n 

Our first test bed is the 30 GHz channel of Planck/LFI. 
We remind here that Planck spins at 1 rpm, has a bore- 
sight angle of 85° and observes the same circle of the sky 
for 1 hour (every hour the spin axis is moved, say along 
the ecliptic, by 2.5'). The mission is planned to remain 
in operation for at least 14 months, resulting in approxi- 
mately two nearly full sky coverages. This corresponds to 
~ 10^ samples for each of the two LFI 30 GHz horns. The 
size of this timestrcam makes Monte Carlo simulations 
quite prohibitive. Therefore, we have decided to use only 
~ 2^'' (~ 2 • lO"^) observations corresponding to ~ 140 
rings (each ring is scanned 60 times) and covering two 
slices, each about 12° degrees wide. This means that we 
process ~ 12 days of observation. For the moment, we 
choose the usual "1//" noise given in Eq. (^). As dis- 
cussed in Sect. 2, we will examine both the a < 1 and 
the a > 1 (i.e., non stationary) case. The values of A and 
fk are chosen accordingly to the Planck specifications: 
A ~ 4 • lO'' ( ^J.K^/s)'^ and fk = 0.1 Hz. The simulations 
are performed by first generating a pure noise timestrcam, 
n, according to Eq. (^ with the power spectrum given in 
Eq. (||). This is accomplished by use of standard EFT 
techniques. The n estimator is trivially computed using 
Eq. (|l4|). The periodograms of the series are then esti- 
mated by binning together a given number of neighbor- 
ing frequencies (in the statistical literature this procedure 
is known as "pooling of the periodogram"). This choice 
has two advantages: (1) it brings down the size of the 
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Fig. 1. LP estimates for a = 0.7 on n (black curve) and Fig. 2. Same as Fig. |l| but for WE. Note the shrinking in 
n (dark blue curve). The bell shaped curves overplotted the horizontal scale. 



are Gaussian fits to the data. 



dataset (from 2^^ proper frequencies to, say, a few thou- 
sands) while preserving most, if not all, of the relevant in- 
formation; (2) it is potentially beneficial from the point of 
view of the efficiency of the estimator (see Robinson, 1995a 
for discussion concerning the LPE case). Simulations show 
that WE and LPE are affected, though not dramatically, 
by the number of frequencies pooled together; in fact, if we 
bin too many frequencies we can even degrade the qual- 
ity of the estimator. We found that a reasonable choice is 
to pool together 64 frequencies. As discussed in Sect. 2, 
LPE only exploits the m lowest frequencies of the pe- 
riodogram. In our simulations, we find that an optimal 
number is m ~ 80.^ While the implementation of LPE 
is straightforward, requiring only a linear regression rou- 
tine, the same is not true for WE. In fact, in the latter case 
we have to minimize Eq. (|l^) over the chosen parameter 
space which is spanned by all physically acceptable values 
of A, fk and a. A minimization routine is then needed. 
A seemingly good candidate is a direction set algorithm 
(see e.g. Powell's method described in Press et al. 1992) 
which is unfortunately too slow for our purposes. A much 
faster choice is to use a variable metric method, such as 
the BFGS algorithm (see again Press et al. 1992). This 
method requires that the gradient of the function of in- 
terest is calculated at an arbitrary point, an information 
readily obtained by differentiating Eq. (^ w.r.t. the pa- 
rameters. 

We have performed 1000 realizations for a = 0.7, 1 and 
1.3, respectively. Our results for LPE and WE are shown 
as histograms in Figs. |l| to ^ and summarized in Table |] 
where mean values and standard deviations (across real- 
izations) are reported. 



In the statistical literature it is occasionally suggested 
to trim (i.e. drop) the lowest (~ 10) frequencies of the peri- 
odogram. We have verified that, in our case, this choice has no 
significant impact on results. 




0.9B 1.00 



1.02 1.04 1.06 



Fig. 3. Same as Fig. hi for a = 1. 



Our first conclusion is that both LPE and WE are not 
significantly degraded when n is considered. In fact, the 
bias is at most 1% for WE and only slightly larger for 
LPE. The results for standard deviations are also compa- 
rable. Indeed the relative performance of LPE increases 
as a grows. In general, WE performs slightly better, as 
expected. However, the simple, robust LPE stands up as 
a very convincing alternative in the absence of precise a 
priori information on the functional form of the spectral 
density. This point will be further discussed below. 

We have also implemented the Moulines and Soulier 
(1999) broad-band log periodogram estimator mentioned 
in Sect. 2. The resulting estimates, however, were ex- 
tremely close to those obtained with LPE and we thus 
decided to omit their presentation for the sake of brevity. 



3.2. Effects on Map Making 

In order to understand how WE and LPE impact on map 
making we evaluated GLS maps out of 50 data streams 
simulated as discussed in Sect. 3.1. We recall here that 



6 P. Natoli et al.: Non- iterative methods to estimate the in-flight 

Table 1. Recovery of relevant parameters for WE and LPE 
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1 
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1. 

n 


3 

n 


LPE 
Qlp 


0.671 
±0.015 


0.657 
±0.016 


1.013 
±0.017 


0.988 
±0.016 


1.361 
±0.021 


1.314 
±0.016 


WE 
aw 


0.7015 
±0.0028 


0.7036 
±0.0029 


1.0042 
±0.0033 


1.0082 
±0.0039 


1.3223 
±0.0054 


1.317 
±0.030 



3 60 - 





Fig. 4. Same as Fig. |[ for a = 1. 



Fig. 6. Same as Fig. |], for a = 1.3. 




estimates of a. So, a fiducial value for fk is used to com- 
pute N. This approach finds its theoretical justification in 
Dahlhaus (1995). He shows that taking 



iV. 



3k OC 



Aj-k)f 



(15) 



provides estimates asymptotically equivalent to exact GLS 
under the only assumption that Eq. (^) holds. Note that 
Eq. ( p^ only depends on a. 

A natural figure of merit to assess the relative perfor- 
mance of LPE and WE is the angular power spectrum of 
the maps. The latter is computed as 



Fig. 5. Same as Fig. |], for a = 1.3. 

the GLS estimate for the map is: 

m= (P'N-ip)-ip'N^M, 

where N is the noise covariance matrix introduced in 
Eq. (0). The Monte Carlo pipeline relies on the follow- 
ing scheme: (1) generation of the noise timestream, n; (2) 
computation of fi; (3) derivation of WE and LPE. The 
maps are computed using the iterative map making algo- 
rithm discussed in Natoli et al. (2001). 

We note that WE return estimates of A, fk and a. 
Only the latter two are required to build N, the scale fac- 
tor being irrelevant. On the other hand, LPE only returns 



2i+l 



(16) 



where the a^m's are coefficients in the spherical harmonic 
expansion of the noise maps. We are not interested here 
in a rigorous power spectrum estimation but, rather, in a 
comparison between LPE and WE. This justifies the use 
of Eq. ( p^ which neglects corrections for, e.g., the sky 
coverage and detailed shape of the surveyed region. The 
results are shown in Fig. ^ (note that the Ces are averaged 
over the simulation index). The curves for WE and LPE 
are virtually indistinguishable. This confirms that only the 
knowledge of a is relevant for GLS map making. Also, the 
small differences between LPE and WE become de facto 
immaterial when map making is at stake. This is an impor- 
tant point: as discussed in Sect. 2, LPE is semiparametric 
and as such less prone to suffer from a misspecified model. 
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Fig. 7. Averaged angular power spectrum of 50 maps dis- Fig. 8. As for Fig. ^but comparing WE (blue curve) with 

played as a function of multipole The light blue line the iterative (hght red) method. The input value for a is 

is for WE and completely hides the curve for LPE. The 0.7 
bottom red line is the difference between the two. The 

spectral index input value is a = 1.2. 



p.d.f. of "medium" observed pixel (660) 



The heuristic motivation behind this result can be eas- 
ily understood if we think of GLS map making as equiva- 
lent to OLS performed on a pre-filtered series. Pre-filtering 
needs not to be based on the exact form of the noise spec- 
tral density. What really matters is to strongly suppress 
the long range correlation and this can be achieved even 
with an imperfect estimate of the noise memory parame- 
ter. 

3.3. Comparison with iterative methods 

The purpose of this Section is to compare the efficiency of 
LPE and WE against iterative noise-signal estimators. It 
is generally claimed (Prunet et al. 2001, Dore et al. 2001) 
that, when reasonable criteria are chosen, convergence on 
the noise estimator is reached fairly quickly, and only a few 
(~ 4 to 5) iterations are sufficient for most applications. 

The following scheme is employed to implement the 
iterative method. We start from the OLS noise estimator, 
n. The iteration is carried out by: (1) using the latter to 
estimate N; (2) performing map making to obtain a new 
solution for the map, rh; (3) estimating the noise again as 
fi = d — Pm. This procedure, basically the same described 
in Prunet et al. (2001) and more recently in Dore et al. 
(2001), has been implemented by using our map making 
tool (Natoli et al. 2001), performing 6 iterations (1 OLS 
plus 5 GLS map making runs) for each map, out of a total 
of 50 Monte Carlo realizations. In the same spirit of Fig. ^ 
we derive CiS for these maps and compare them with the 
results based on WE (this time for an input value a — 0.7). 
Results are shown in Fig. ^: the similarity between the two 
curves is striking. 

To further compare iterative and non iterative meth- 
ods we have computed histograms showing the frequency 
of occurrence of temperature values for selected pixels of 




0.0 
ink 



Fig. 9. Empirical distribution of the "medium" hit pix- 
els (660 hits). WE (triangles) and the iterative scheme 
(squares). 



the maps. We have chosen a "medium" observed pixel, 
which is hit 660 times in the Planck scan described above 
and a "low" observed pixel, hit 240 times. The histograms 
are shown in Figs. ^ and ^ respectively. Note that all 
the pixels with the given hit number contribute to the 
histograms which do not show outliers. 

It should be clear at this point that estimating iter- 
atively the noise properties yields no obvious advantage 
when compared to the simpler WE and LPE described 
in this paper. Furthermore, as stressed in Sect. 1, the 
iterative scheme described above needs a sound statisti- 
cal foundation in the presence of long range dependence. 
Moreover, it is important to realize that such a scheme is 
severely more costly than LPE and WE. In fact, in order 
to implement the latter, only a one step OLS map mak- 
ing is needed, plus one periodogram estimation over the 
time line and the final GLS one shot map making (the ac- 
tual computational effort to derive WE and LPE from the 
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Fig. 10. Same as Fig. ^but for pixels observed 240 times. Fig. 11. Same as Figs. I and 0, but for a BOOMERanG- 

like simulation. The histogram are "observed" distribu- 
tions of GLS mapmaking temperatures for all pixel hit 23 

periodogram is extremely small). On the other hand, the times in the maps. Noise has been estimated using WE 

iterative method requires, other than the OLS map mak- (black curve) and LPE (dark blue). 

ing, a periodogram estimation and a GLS map making at 

every noise iteration. ^^g^ 



3.4. Application to BOOIVIERanG 

All simulations shown above were performed for a sub- 
set of the Planck timestream. Furthermore, the noise 
model assumed as input is the customary "1//" noise (see 
Eq. (^). In this Section we relax both assumptions and 
apply our techniques to BOOMERanG. 

BOOMERanG is a balloon borne experiment which 
scans at constant elevation by allowing the gondola to 
slew at 1 or 2 degrees per second (d.p.s.), while the detec- 
tor sampling rate is 60 Hz (de Bernardis et al. 1999). To 
make our simulations realistic we take a chunk (2^^ sam- 
ples out of total of about 3 • 10^) of the 1 d.p.s. part of the 
scan performed by the B150A channel during the flight 
of 1998. Contrariously to Planck, this scan is aimed at 
a small patch. The central region comprises a few tens 
of square degrees and its coverage (i.e. number of hits 
per pixel) is far more uniform. The noise properties of 
the BOOMERanG detectors are not well described by the 
model given in Eq. (|^) because the timestream is altered 
to deconvolve bolometric filters (both low pass and high 
pass). To take these effects into account we propose to use: 




-500 

Pixel temperature i/^K) 



P(/) = Asin^(|^ 



1 



(17) 



with fiducial values A ~ 2 • 10"* ( fiKy/s)^, fo ~ 10"^ Hz, 
/i ~ 7 Hz, fk = 0.15 Hz, A' - 10« and a = 1.5. We choose 
fo fi, since the exponential term arises after deconvon- 
volving the detector high frequency response. To avoid 
excessive contamination from this high frequency term we 
sharply cut the spectrum at 20 Hz (as in BOOMERanG: 
see Hivon et al. 2001). 

We thus use Eq. (||) to generate Gaussian noise assum- 
ing the power spectrum given in Eq. (p7|). Subsequently, 



Fig. 12. Same as Fig. ^ but for a pixel hit 100 times. 
We compare here WE (dark blue curve) with the iterative 
method (black). The light blue curve is computed using 
the correct noise power spectrum (Eq. (|l7|)) and, hence, 
provides the exact GLS estimate. 



WE and LPE are applied to the data. WE is fully para- 
metric and requires that a model for the spectral density is 
specified. Rather than supply the model given by Eq. (|l7|), 
used to generate the data and hence "exact" , we base the 
estimations on Eq. (^, i.e. on standard "1//" noise. This 
results in a deliberate model misspecification, as we are 
providing WE with an incorrect model. This choice has 
the obvious purpose to test the robustness of the esti- 
mate. Strictly speaking, for a < 2 Eq. ( [l7| ) does no longer 
display long memory behavior, and in this sense we are 
misleading LPE as well. However, long memory is still a 
good approximation if fo ^ fk- 

Our results are shown in Fig. |ll] where we give, for 
the most frequently observed pixel in our scan, temper- 



p. Natoli et al.: Non- iterative methods to estimate the in-flight 



9 



ature histograms analogous to the ones given in Figs. ^ 
and The similarity between the two curves is striking. 
Very much the same conclusions as in the previous Section 
can be drawn: the difference between WE and LPE, and 
between WE and the iterative method, is by all means 
negligible as shown in Fig. In the same Figure we plot 
the exact (unfeasible) GLS estimate. Again, no significant 
difference emerges. 

4. Summary and conclusions 

This paper has addressed estimation methods for the 
properties of CMB receiver noise which is well approxi- 
mated by a long memory process. The most recent litera- 
ture proves that optimal statistical inference is deeply af- 
fected by such characteristics. Under these circumstances 
results by Dahlhaus (1995) entail that map making proce- 
dures can be performed with a single iteration, provided 
that the slope at the origin of the log spectral density (i.e. 
a) is known; such procedures are asymptotically equiv- 
alent to maximum likelihood, GLS-type estimates. Our 
prime effort, hence, has been to address the efficient es- 
timation of the memory parameter a. We have first pre- 
sented the fully parametric WE which allows for the full 
correlation structure of the noise process to be recovered. 
WE, however, requires that the functional form of the 
noise spectral density is known. To relax this very restric- 
tive assumption we also focus on the semiparametric case, 
where it is only assumed that a spectral singularity be 
present at the origin. Both the parametric and the semi- 
parametric procedures are known to enjoy some nice ro- 
bustness properties in non-standard situations, including 
non-Gaussian and non-stationary circumstances. 

We assess the performance of these statistical proce- 
dures by implementing them on simulated Planck and 
BOOMERanG datasets. We note first that for practical 
implementation we need to filter out the signal component 
by running a preliminary OLS regression. We have shown 
by direct Monte Carlo comparison that this induces no 



correct evaluation of the dependence structure of receiver 
noise may in principle improve the determination of con- 
fidence intervals for the angular power spectrum and its 
governing parameters, either by resampling methods or 
by asymptotic approximations. These issues are currently 
being investigated. 
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